Reduced equations of motion for quantum systems driven by diffusive Markov 

processes 



Mohan Sarovar^ and Matthew D. Grace^ 

^Department of Scalable and Secure Systems Research, 
Sandia National Laboratories, Livermore, CA 94550 

The expansion of a stochastic LiouviUe equation for the coupled evolution of a quantum system 
and an Ornstein-Uhlenbeck process into a hierarchy of coupled differential equations is a useful 
technique that simplifies the simulation of stochastically-driven quantum systems. We expand the 
applicability of this technique by completely characterizing the class of difi^usive Markov processes 
for which a useful hierarchy of equations can be derived. The expansion of this technique enables the 
examination of quantum systems driven by non-Gaussian stochastic processes with bounded range. 
We present an application of this extended technique by simulating Stark-tuned Forster resonance 
transfer in Rydberg atoms with non-perturbative position fluctuations. 



Describing the dynamics of a quantum system coupled 
to uncontrolled degrees of freedom has been an impor- 
tant problem since the inception of quantum mechanics, 
e.g., [T]. The accurate description of such open quantum 
systems is particularly vital for the design of quantum 
technologies, such as quantum computers. Several ap- 
proximate and exact methods exist for describing the 
dynamics of open quantum systems, including master 
equations, surrogate Hamiltonians, and Monte-Carlo nu- 
merical simulations. 

In this work, we examine quantum systems driven lin- 
early by classical fluctuations. That is, the Hamiltonian 
for the system is described by Hn{t) = HQ-\-Q.{t)V , where 
U,{t) is a time-dependent stochastic variable. This is a 
sub-class of a more general open quantum system where 
ri(t) is an operator in the Hilbert space of an uncon- 
trolled environment. The replacement of the operator 
with a scalar variable is an approximation that is valid 
in certain limits (e.g., in the high temperature limit of 
a bosonic environment). Scalar fluctuations can also de- 
scribe noise in a quantum system that is controlled by an 
effectively classical quantity, such as a gate voltage. 

For such stochastic evolution, the dynamics of the sys- 
tem under a given trace of the noise is dictated by the 
von-Neumann (Schrodinger) equation: 

^mmt)}) = -\[Hn{t),m{m})]- w 

Here, the notation p{t\{Q.{t)}) is used to explicitly in- 
dicate that this density matrix is conditioned on a par- 
ticular realization of the noise process. This equation is 
formally solved to yield p{t\{n{t)}) = U{t,Q)fPU{t,Q)\ 
for initial state where U{t^G) = Te^s-'n^"^*) with 
T being the time-ordering operator. We are often more 
interested in the unconditioned evolution of the system 
state, after the fluctuating quantity has been averaged: 



p{t) = (/5(t|{f](t)})){o(t)}, 



(2) 



where the angled brackets denote an expectation value 
over the stochastic process up to time t. We will re- 
fer to differential equations describing the evolution of 



this averaged quantity as reduced equations of motion. 
In this work we derive reduced equations of motion for a 
quantum system coupled to a wide family of stochastic 
processes. 

We focus on stochastic processes that are time- 
continuous, time-homogenous, and Markov, which means 
that the evolution of the conditional probability distribu- 
tion for the process evolves as 



d_ 
dt 



p{n,t\n' ,t') ^TnP{n,t\ri' ,t'), 



(3) 



where P{il,t\il' ,t') is the probability that the stochastic 
process takes the value at time t given that it took 
the value fi' at time t' {t' < t). is the forward gen- 
erator of the process and is a differential operator only 
involving derivatives with respect to fl. The generator 
of evolution is in general a complex quantity that results 
from a Kramers-Moyal expansion of the (classical) mas- 
ter equation for the probability distribution of the process 
J2|. Below, we will restrict ourselves to certain forms of 
this generator, but for now only the Markov assumption 
will be used. 

In order to derive an expression related to the reduced 
equations of motion in a simple manner, consider the 
quantity p(t,ft) = pfi{t)P{fl,t), which is the joint distri- 
bution of the quantum-mechanical coordinates and the 
stochastic variable at time t. pn{t) = p{t\n,{t)) is the 
density matrix conditioned on the value of the stochas- 
tic process at time t (in contrast to p{t\{n{t)}), which 
is the density matrix conditioned on an entire history of 
the stochastic process) . Evaluating the time derivative of 
this joint distribution results in the stochastic LiouviUe 
equation, flrst derived by Kubo [5]: 



dt 



':Q.{t)V^ 



p{t,n) + Tnp{t,^l^) 



where B = [A,B]. Here, we have used the fact that 
because the stochastic process is Markov, the time evo- 
lution of P(ri, t) follows the same law as the conditional 
distribution in Eq. ([s]) and is generated by Tq. This 
simple linear form for the evolution of the joint distri- 
bution is possible because both the conditional quantum 
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density matrix and the probability distribution for the 
stochastic variable evolve linearly and in a Markov fash- 
ion. Because /5(i, fi) is a joint distribution, its marginal 
over the stochastic variable yields the average density 
matrix, i.e., p{t) — J^^ d^lp{t, fl), where T) is the range of 
the stochastic variable [6]. 

Markov diffusion processes. A wide a class of phys- 
ical processes can be approximated by truncating the 
Kramers-Moyal expansion for the generator Tq of the 
Markov process at the second term P]- This results in 
the Fokker-Planck equation for the probability distribu- 
tion, with the generator being the differential operator 
a: r^P ^ = -iiA{n) + \ipB{n), where A(n) 
and B{Vt) are real differentiable coefficients, with the re- 
striction that B > 0. 

For such diffusion processes, it is common to define 
the backward, or adjoint, generator for the process: 

= + h^^^)wn- ^ Markov process, 

Eq. ([3]) is often called the forward Kolmogorov equation 

and ^P{n,t\n',t') = %<:i>P{n,t\n',t') is the backward 

Kolmogorov equation. Let /„(0) and be eigen- 

functions of the forward and backward generators, re- 
spectively, i.e., ^nfn{^) = -Kfn{^) and $0 /„(r2) = 
— with A„ > 2J. These two eigenfunctions 
are related through the stationary distribution of the pro- 
cess, explicitly /„(f2) = P°(f2) / „(ri), where P° is the 
stationary distribution of the process fi 3 . Furthermore, 
the eigenfunctions of the backward generator form a com- 
plete, orthogonal system (on the range V) under the mea- 
sure induced by the stationary distribution of the process: 
dVlP^lVl) f „(ri) / m(ri) = Smnj which in turn implies 
the following about the eigenfunctions of the forward gen- 
erator |3]: dQiP^ (fl))-^ Ui^) = S,nn 

Diffusive hierarchical equations of motion. Us- 



ing the completeness and orthogonality of the backward 
and forward generator eigenfunctions 0, we can ex- 
pand t) in terms of these eigenfunctions, and conse- 
quently, p{t,n) = pn{t)Pin,t) = Y.n=Qpn{t)fnm Here, 
Pn{t) are unnormalized auxiliary density matrices, which 
do not have the interpretation of being conditioned den- 
sity matrices of the quantum system. Then, expressing 
the stochastic Liouville equation Eq. Q in terms of this 
eigenfunction expansion yields 

00 00 . 

g-^Pn{t)fnm = -'-H^p^{t)f,,{n) 

-'-V'pnit) [n{t)fnm - x„pn{t)fnin). (5) 

To simplify further, we make the assumption that the 
eigenfunctions of the backward Markov generator are 
polynomials in il, i.e., is an n*^ order polyno- 

mial in 17. With this assumption, we utilize the following 
theorem, typically attributed to Favard [7j, which states 
an important property of orthogonal polynomials: 

Theorem (Favard): Let {pn{x)},n > be a system 
of polynomials. This system satisfies a three-term recur- 
rence relation Pn+i{x) = {A^x + Bn)pn{x) - CnPn-l{x), 
(with p-i{x) = 0) if and only if it is a system of or- 
thogonal polynomials. Here An ^ 0, Bn and C„ ^ are 
fi-independent (real) recurrence coefficients. 

Favard's theorem implies that the orthogonal polyno- 
mial eigenfunctions of the backward Markov generator 
(and hence, those of the forward Markov generator) sat- 
isfy the three-term recurrence relations: 

nf^m = ^fn-im - ^fnm + -^fn+im. (6) 

Using this recurrence relation, the eigenfunction expan- 
sion of the stochastic Liouville equation becomes: 



J2 g-/nm„m ^ Y 



Tl = 



n=0 



Pn{t)fnm- -V'pnit) 



^fn+m + ^fn-m - ^frM 
-^n -^n -^n 



(7) 



Multiplying both sides of this equation by f„i{Vl)/P^{Vl) and integrating over SI results in an equation for each m 
because of the orthogonality of the functions fn{^)- These equations form a hierarchy of coupled (operator) differential 
equations, which we refer to as diffusive hierarchical equations of motion (DHEOM) j8j: 



d_ 
dt 

dt 



Po{t) = - 

Pn{t) = - 



Bo - -^y 



+ Ao 
+ A„ 



Po(t) - 

Pn{t) 



i Ci 



V^-piit), 



hAi 

hAn+i nAn^i 



71 > 0. 



(8) 



These equations are useful because they describe the evo- 
lution of the system density matrix under the infiuence of 
stochastic noise, but without explicit reference to noise 
variables. In addition, each member of the hierarchy only 



couples to its two neighbors, i.e., n couples to ?i + 1 and 
n—l, making their integration easy. We also need to spec- 
ify the initial conditions and the prescription for calcu- 
lating the system density matrix from the solution of the 
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hierarchy of coupled equations. A physically motivated 
initial state for the quantity Jl) is p{0, il) = fPP'^{fl), 
where P*^ is the stationary (equilibrium) distribution of 
the noise process and pP is any system density matrix. 
Evolution from this initial state describes the response 
of a quantum system to fluctuations around the bath 
equilibrium. Now, ^qP'^ = by definition, and there- 
fore, po{0) = /5°, and Pn(0) = V n 7^ 0. At any time, 
the system density matrix is defined as the integral of the 
quantity p{t, fi) over fi: 

p p oc 

p{t)^ / dnp{t,a)^ / dl7 Vp„(t) cxpo(t). 
Jv Jv 



The proportionality is a result of the following property 
of the eigenfunctions: 



dQ /„(r2) cx 



V 



4- 



where this proportionality follows from the fact that / ( 
is a zeroth-order polynomial in £7. As we shall see below, 
it is always possible to choose / = 1, and therefore, the 



proportionality above can be converted to an equality. 
Hence, the auxiliary density matrix po{t) always keeps 
track of the averaged system density matrix p{t). 

To summarize, we have shown that the dynamics of 
a quantum system that is linearly driven by a diffusion 
process with a Markov generator possessing polynomial 
eigenfunctions can be described by a hierarchy of cou- 
pled dynamical equations. The solution to these equa- 
tions will reproduce the reduced density matrix of the 
quantum system at any time, with no approximations. 
However, the above hierarchy of equations is infinite, and 
therefore, we require some truncation strategy in order 
to solve them. If the equations can be truncated at some 
n ^ N at the expense of bounded error, they can be nu- 
merically solved (or analytically solved via a partial frac- 
tion expansion if Hq is time- independent j9j). A general 
truncation strategy exists if the following conditions hold: 



(1) for large enough n = N, \\n\ ^ \\Hq 



A. 



and (2) 



where 



1 2 is the induced 2- 



norm and /„ € oj{gn) is indicates that /„ dominates (?„ 
asymptotically (in n) . In appendix |b] we develop such 
a general truncation strategy that results in the follow- 
ing terminator equation for the hierarchy at a level N 
satisfying these conditions: 



d_ 
dt 



PN{t) = - 



rrX ^Wt/> 



+ A 



N 



PN{t) - 



1 



c 



N+1 



-V'V'' -pN{t) - 



1 



HAn-i 



V^pN-iit). (9) 



We now turn to characterizing the diffusive stochastic 
processes with generators that have polynomial eigen- 
functions. A quantum system driven by each one of 
these processes will have a dynamical description in terms 
of the DHEOM derived above. The combination of 
Bochner's theorem [TU], which is stated explicitly in ap- 
pendix |Xj and Favard's theorem implies that there are 
only three diffusive processes with backward and for- 
ward generators that have orthogonal polynomial eigen- 
functions: (a) the Ornstein-Uhlenbeck (OU) process, de- 
fined on (— C)o,C)o) (b) the square-root process, defined 
on (r2niin,cx)) and (c) the Jacobi process, defined on 
(ilmi„, ilniax)- In table [l] of appendix |Aj we list these dif- 
fusive processes with their orthogonal polynomial eigen- 
functions and other relevant properties. In appendix [C] 
we also explicitly write the DHEOM that describe the 
dynamics of a quantum system linearly driven by each of 
the above three stochastic processes. The Jacobi process 
is particularly efficient to simulate because the eigenval- 
ues A„ of its generator scale quadratically with n, making 
truncation possible at a smaller depth than for the other 
two processes, where A„ scales linearly with n. We note 
that the DHEOM for the OU process was previously de- 
rived by Kubo and Tanimura [51 [U] , but this more general 
framework resulting in DHEOMs for all major diffusive 
Markov processes was heretofore unknown. For reasons 



explained in appendix [Cj the DHEOM method is limited 
to simulating square-root processes with mean reversion 
rate 7 > 1, whereas any 7 is valid for the other pro- 
cesses. Finally, for completeness appendix [D] presents 
the DHEOM for the propagator for the dynamical map 
as opposed to the density matrix. 

Application and Discussion. The primary advan- 
tage of our general formulation of DHEOM is that it 
allows for the efficient simulation of quantum systems 
driven by noise sources with bounded range and non- 
Gaussian amplitude distribution. As an application, we 
consider energy transfer between two systems mediated 
by a dipole-dipole interaction. This interaction is an ap- 
proximation to the Coulomb coupling of electric charge 
distributions and is the basis of several common exper- 
imental techniques, such as Forster resonance energy 
transfer (FRET) spectroscopy [TT] and Stark-tuned reso- 
nant transfer in Rydberg atoms [12l HI] ; we focus on this 
latter system as our example. 

The electric dipole-dipole coupling energy scales as 
J (X l/i?"^, where R is the distance between the coupled 
molecules or atoms. Consider the case where this inter- 
particle spacing is varying in time as a result of ther- 
mal, electromagnetic, or vibrational fluctuations. If the 
fluctuations in R are treated perturbatively the dipole- 
dipole coupling energy can be expanded linearly around 
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FIG. 1. Average population in the excited state of atom 2 
after an interaction time of T = l^s as a function of the Stark 
detuning for different models of motional noise. The black 
(solid) curve is for coherent evolution with no fluctuation in 
inter-atom distance. The parameters used are: 7 — 1.5, /i = 1 
for all noise sources, and = 0.3 for the OU process. The 
square root process is defined in the semi-interval [0, 00) with 
ci — 1, and the Jacobi process in the interval (1/8, 8) with c = 
1. See appendix [A] for the interpretation of these parameters. 



the mean separation, i?o : J ^ l/Ro — {3/Rl)SR{t). How- 
ever, in cases where such a perturbative treatment is 
not accurate we must consider the energy J as a time- 
dependent fluctuating quantity in a range [Jmin, >/max]- 
To explore the consequences of such a non-perturbative 
treatment, consider the Hamiltonian of two interacting 
Rydberg atoms: H = —eial — e2crl + J{t){(T]_a^+ahcr'^), 
where we restrict the description of the atoms to the 
two energy levels relevant to the energy transfer, and 
therefore use Pauli matrices to describe them as effective 
two- level systems in the basis |r!_) , |r!|_) where r^. are 
Rydberg states for atom i with the energy gap e^, i.e., 
CT* |r!|_) — ±ei |r!|_). Typically, these levels are chosen 
so that when an electric field tunes atomic energies by 
the Stark effect they satisfy ei = €2 (Forster resonance 
condition) at some critical field value. For example, 
in Ref. m, K) = \37Py2) ,\rl) = |375i/2) , |r^) = 
|38S'i/2) , |r^) = I37P3/2), where these are fine states of 
Rb atoms. We ignore relaxation of these states because 
for sufficiently cooled Rydberg atoms it can be neglected 
on the timescales we are considering [2] . We also assume 
that the relative orientation of the atomic dipoles remains 
constant. The initial state is the excited state of atom 
1: p° = |rjj_, r^) (rj'j_, r?. |, and the dipole-dipole coupling 
is assumed to fluctuate as J{t) = Jofi(i), where il{t) is 
one of the three diffusive noise processes described above. 
The quantity being observed (e.g., by selective field ion- 
ization [13] ) is the population in the excited state of atom 
2 after an interaction time T. Although position fluctu- 
ations were reported to be minor in the experiment of 
Ref. |13| . we use physical parameters compatible with 
this experiment for concreteness: Jo/fi = 0.5MHz, and 



T = 1/is. Fig. [T] shows the average population of |r^) at 
T for different values of the Stark detuning A = 62 — ei . 
We expect that this transferred population will drop off 
as |A| increases (as the atoms become less resonant), but 
the figure shows that the behavior of this drop-off de- 
pends heavily on the exact nature of the fluctuations in 
R, although all processes have the same mean fJ. = 1 and 
mean reversion rate 7 = 1.5. The OU and square-root 
processes predict average populations that have signifi- 
cant oscillations as a function of the detuning, similar to 
the completely coherent (noiseless) case. In contrast, the 
Jacobi process driven evolution predicts damped oscilla- 
tions and smaller transfer populations. We predict that 
for experiments performed in this regime, where position 
fluctuations are relevant, the Jacobi process average will 
be more accurate (and they are certainly more consistent 
with what is experimentally observed in Ref. [13 ). This 
is because the Jacobi process is a more accurate descrip- 
tion than the OU and square-root processes, which can 
only be approximations in this scenario; realistic trapping 
conditions impose strict lower and upper bounds on J re- 
sulting from the fact that the atoms can only drift within 
the trap volumes. The OU and square-root simulations 
predict large transfer populations because rare, unphysi- 
cal, large magnitude couplings contribute to the averages 
in these cases while they do not to the Jacobi process av- 
erage. These results emphasize the advantage of being 
able to efficiently average over non-Gaussian noise pro- 
cesses with bounded range. To underscore the numerical 
efficiency of the DHEOM simulations, we note that the 
calculations producing Fig. [T] were performed 10 times 
faster using the present theory than conventional Monte 
Carlo sampling of the noise processes |15] . Furthermore, 
because of the well-defined truncation strategy, conver- 
gence issues and noise are not a concern for DHEOM 
simulations as in the Monte Carlo approach. 

Summary. We have completely characterized the 
class of diffusively driven quantum systems for which a 
useful hierarchy of reduced equations of motion can be 
derived, and explicitly derived these DHEOM. We ap- 
plied the technique to examine position fluctuation de- 
pendence of Stark tuned Forster resonance between Ryd- 
berg atoms. We expect that the general techniques devel- 
oped in this work will be useful for widening the range of 
open quantum systems that can be simulated efficiently. 
Use of the DHEOM for Jacobi processes will be partic- 
ularly useful for modeling quantum systems subject to 
fluctuations with bounded range, e.g., classical noise in 
external control flelds [T^. 
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Appendix A: Bochner's theorem and polynomial eigenfunctions 

In the main text, we claim that there are three stochastic processes that admit an eigenfunction of their Markov 
generator that is an orthogonal polynomial. This result follows from Bochner's theorem [TU]: 
Theorem (Bochner): Consider the eigenvalue problem 

(P d 

/(2:)^y(a;) + 9ix)-^yix) + Hx)y{x) = Xy{x) (Al) 

where f,g,h are fixed polynomials and A is independent of x. There are only five classes of polynomial solutions to 
this problem, where yn{x) is a polynomial of degree n, and they are: 

1- yn{x) = He„(a;),n £ Nq is a Hermite polynomial, defined on the infinite interval x G (— cx3,oo). 

2. yn{x) ~ Li^\x),n G Nq is a generalized Laguerre polynomial, defined on the semi-infinite interval x G [a, oo) 
or a; G (— oo, a] for a G M. 

3- ynix) = ji"'^\x), n G No is a Jacobi polynomial, defined on the interval x G (a, b) for a, 6 G M. 

4- yn{x) is in the set {x"} for n G Nq. 

5. ynix) — Bn{x) is a Bessel polynomial, defined on the infinite interval x G (—00,00). 



Note that the backward Kolmogorov equation is an example of the differential equation considered in Eq. ( Al ), with 
the additional restriction that f{x) > 0. Therefore, all polynomial eigenfunctions of the backwards Kolmogorov 
generator fall within the Bochner classes identified above. However, since we also require that our eigenfunctions 
satisfy a three-term recurrence relation Favard's theorem further restricts us to the subset of the Bochner classes that 
are orthogonal polynomials. The Hermite, Laguerre, and Jacobi polynomials are all classical orthogonal polynomials 
and therefore obviously orthogonal. The polynomials sequence in class 4 above is not an orthogonal sequence, and 
similarly, the Bessel polynomials are not an orthogonal set of polynomials for any positive-definite measure. The Bessel 
polynomials are orthogonal under a quasi-definite measure and hence do satisfy a three-term recurrence relation [10] . 
However the positive-definiteness of the measure is an important feature in the following and therefore we will not 
consider the Bessel polynomials. Therefore, the combination Favard's and Bochner's theorems imply that there are 
only three diffusive processes with backward and forward generators that have orthogonal polynomial eigenfunctions: 
(a) the Ornstein-Uhlenbeck process, (b) the square-root (or Cox-IngersoU-Ross) process, and (c) the Jacobi process. 
This result has been derived by a number of authors in the past, in different mathematical contexts, including in Refs. 



Appendix B: Truncation of DHEOM 



The conditions for a general truncation strategy to exist for the DHEOM are: (1) for large enough n = N, 

\Xn\ 3> \\Hq — ^V^\\2, and (2) ^ G w (^x^) ; where || • II2 is the operator 2-norm and /„ G uj{gn) is indicates that 

/„ dominates gn asymptotically (in n). To develop the truncation strategy under these conditions we follow Ref. 
and consider the formal solution to the n = + 1 ^ 1 member of the hierarchy: 



PN+l{t) 



dr exp 



B 



N+l 



A 



N+l 



N+l 



(t-r) 



The first condition above implies that by choice of N we have |AAr+i| ^ 





Cn+2 




Ho- 


^N + l T 



Pn+2{t) ■ 



1 



(Bl) 



jV\\, and we can replace pi{T) 



with pi{t) in the integrand because the exponential term decays much faster than the rate at which pi{T) changes. 
This approximation enables the integral to be computed and results in the solution 



PN+lit) 



hx 



N+l 



Cn+2 - , s 1 _ , . 

PN+2{t) + -r-pN[t) 

An+2 An 
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TABLE I. Summary of diffusion processes with backward generators possessing polynomial eigenf unctions. I(a,b) is the indicator 
function for the interval (a, 6), r[-] is the Gamma function, and B[-] is the Beta function. All processes have exponentially 
decaying correlation functions, proportional to e"''*. Note that in order to get convenient recurrence relations, we choose 
unconventional normalizations for the orthogonal polynomial eigenfunctions. 



Ornstein-Uhlenbeck process 



Range of stochastic process 

■Dou : — oo < f2 < oo 
Generator parameters 

Eigenfunctions of backward generator 

Scaled Hermite polynomials 

Eigenvalues of backward ^nerator 



(-1)" 

Z^m=0 m\2^(n—2m 



)t(v1*("-'^))' 



$H"^„(fi) = -n7 0„(n) 
Stationary distribution 

Gaussian distribution 



A„ = 717 



pOU 



\/ 7rcr^/7 



Square-root (or Cox-Ingersoll-Ross) process 



7 > 0, /u e •PsR 



Range of stochastic process 

©SR : < n < 00, ci>0 
Generator parameters 

ASR(n) = -7(n-/.), 

B^^{^) = cifi + Co 
Eigenfunctions of backward generator 

Scaled generalized Laguerre polynomials 

Vn(fi) = ^^-I'i ' [^(^ + ^)] = ^\ Z)m=0 (n-m) ' rnT^^ 



where a ; 



27 



1 



Eigenvalues of backward generator 

Stationary distribution 

Gamma distribution Po^(f2) = 



= n7 
r[a+i] 



Jacobi process 

Range of stochastic process 

©J : wi < n < u;2 
Generator parameters 

= -7(0-^), 7>o,^e©j 
^■'(n) = -c(n-tJi)(n-co'2), 00 

Eigenfunctions of backward generator 
Scaled Jacobi polynomials 

Y^n /n \ r[a + /3+n + m + l] / H — \ 

^ n! Z-/m3=0 \m/ r[a + m+l] I aJ2— cJi ) 

C C>J2— Wl 



where a ; 



r[a + m+ 

^^^^^ - 1, and P 

c UJ2—^1 

Eigenvalues of backward generator 

= (-n7 - |cn(n - l))tn(n); 
Stationary distribution 



Beta distribution Pq (f2) = 



An = W7 + |cn(n - 1) 



t/c-1 -'(wi,W2) 
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C 1 

Now the second condition above implies that for large enough TV, ^ < ^ and we can ignore the contribution 
from the first term in the square bracket, resulting in: pAr+i(t) h\i^\AN PN{t)- Using this approximation, we 
can terminate the hierarchy at level N and formulate the terminator equation as: 



d_ 
di 



PN{t) 



B 



N 



PN{t) 



1 



c 



N+1 



h? Xn+iAn+iAm 



1 



hA 



N-l 



-V'pN-iit) (B2) 



Appendix C: Explicit forms of the DHEOMs 

In this Appendix we explicitly write DHEOM for quantum systems driven by the three diffusive Markov processes 
that permit such descriptions. We have confirmed that these DHEOM descriptions of dynamical evolution agree with 
averaged Monte-Carlo simulations of the same noise processes. 

1. Ornstein-Uhlenbeck process 

For the case where fl(t) is an Ornstein-Uhlenbeck process, Kubo and Tanimura [5, ^ have demonstrated how to 
construct hierarchical equations of motion for the reduced density matrix. In fact, the generalization described in 
this paper was originally motivated by this work. In this section we restate the Kubo- Tanimura DHEOM for the OU 
process in the present notation. 

The joint distribution is first expanded in terms of the eigenfunctions of the Markov generator for the Ornstein- 
Uhlenbeck process: 



p{t, fl) = pn{t)Pi^, = Pnit)M^) 



(CI) 



n=0 



(j)n{^) = -P(P^(^^) <!> n{^) are eigenfunctions of the forward generator for the OU process. These quantities are defined 
explicitly in Table |l] Since (/)„(r2) are proportional to Hermite polynomials a three term recurrence relation they follow 



(C2) 



where Table |T] should be referred to for the interpretation of the parameters in this expression. Given this recurrence 
relation, the resulting DHEOM is: 



|p-oW = {Ho^ + PVX) p,{t) -'-M v^Mt), 



d_ 

Ft 



Pn{t) 



■ nj 



2 ' / 2 

Pn{t)^'-J'^V''pn+i{t)^jnJ'^V''-pn^i{t), n>0 (C3) 



The terminator equation is explicitly 
d 



dt 



PN{t) 



P^W - ^ (^) V-V-pMt) ~ 'j^Nj^V-pN-iit) (C4) 



2. Square-root process 
We will expand the joint distribution for the quantum and stochastic variables as 

oo 

p(t,l]) = ^p„(i)x„(r!) 



(C5) 
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where Xni^) = P|^(^)V«(f^) is the expression for the eigenfunctions of the forward generator for the square-root 
process in terms of the quantities in Table |l] The three-term recurrence formula satisfied by the basis functions Xni^) 
(which follows from a three-term recurrence relation followed by the Laguerre polynomials) that we require is: 



27 n 

Ci 



Co 



27 



27 Cl 



(C6) 



where a = ^ (^fj, + ~ I, and j,fi,CQ, and ci are the defining parameters for the process (see Table |lj). Using this 
recurrence relation, the resulting hierarchy of coupled equations of motion are: 



= 4'"" 



Cl f , ^ Co 

— a + 1 

27 Cl 



V-]po{t)+'-^{a + l)V-p,{t), 



d_ 
di 



Pn(t) = - 



|1(..2..1)-^ 



- 717 



The terminator equation is explicitly, 
d 



dt 



PN{t) 



m 



^(a-f2iV+l)-^ 

27 Cl 



Pn(t) 



i Cl (a + n + l)_^^^_ . . i ci 



n>0 



PNit) 



(C7) 



^ (a + iv + Dv^v^Mt) + i^N'y'PN-iit) 



(C8) 



Note that this choice of recurrence relation leads to the coefficient of scaling as — > n/7 for large n. At the 
same time, the eigenvalue of the generator scales as — )■ 771. This means that the first condition for the existence of 
a truncation strategy, that |A^r| ^ \\Hq — ^V^\\2 can only hold true for 7 > 1. Therefore we are only guaranteed 
that a valid truncation exists when 7 > 1 for the square-root process. It is possible that this condition can be lifted 
by using an alternate recurrence relation for the generalized Laguerre polynomials but we have not been able to find 
such a recurrence relation at this time, and therefore are limited to simulating square-root processes with 7 > 1 using 
the above DHEOM. 



3. Jacobi process 



If the stochastic driving of the quantum system is by a Jacobi process, we expand the joint distribution in terms of 
the eigenfunctions of its generator 



(C9) 



n=0 



where Cn(^^) = Poi^) C ni^) in terms of the quantities defined in Table |lj The three-term recurrence formula we will 
use, which follows from a three-term recurrence formula of Jacobi polynomials, is: 



nun) = 



Aio f (a + n){/3 + n) 
nr]n{rin + 1) 



Cn-l(f^) 



Au 



UJ2 



h " 

VniVn + 2) 

Auj / {n + iy{r]n-n+l) 
~Y \ (ry„ + l)(r7„ + 2) 



Cn(f^) 



(CIO) 
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where a and /3 are as defined in Table [T] and Ao; 
derive the following DHEOM: 



^Po(i) 



UJ2 



UJ2 



Aid 



Voivo + 2) 



CJ2 ^ i^i , and 77„ = a + /3 + 2n. Using this recurrence relation, we 



1 



?7„(r7„ + 2) 
{a + n + l){l3 + n+l) 
/i""' V ("+ l)'7n+l(??n+l + 1) 
_ / n^iVn-i -n + 2) 

ti V(?7«-i + l)(?7n-i+2) 
The terminator equation for this DHEOM is explicitly, 



1 



V 

V 



_ I Au /(a + lK^+l) 



{nj + -cn[n — 1)) 



71 > 



Pnit) 



g-fPN{t) = 



UJ2 



Auj 



1 



+ (iV7+ -c7V(7V- 1)) 



-Auj 



-Au? 



IN) 



VNiVN + 2) 
{a + N + l){P + N + l)(r;jv - iV + 1) 

VN+iiVN + + + 2) 

-iV + 2) 



PN{t) 



(?7a'-i + 1)('7w-i+2) 



(Cll) 



(C12) 



It is interesting to note that the growth of the A„ with n for the Jacobi process is quadratic (as opposed to linear 
for the OU and square-root processes). This implies that the termination conditions for the DHEOM are met much 
more quickly for the Jacobi process, and hence it is more efficient to simulate than the other processes. 



Appendix D: DHEOM for the integrated map 



In the main text, we derived a general DHEOM for the density matrix describing a quantum system driven by a 
diffusive Markov process. In solving for the dynamics of a quantum systems, it is often useful to work at the level 
of an integrated (completely positive) map, or propagator, rather than a density matrix. The two are related simply 
as p{t) — £{t)p'^, where £{t) is the integrated map that propagates any initial state p° to p{t). If one chooses to 
vectorize the n x n density matrix into an x 1 vector by stacking the columns: ~f} = vec(/5), then £(t) is a matrix 
of size x n^. For unitary evolution, where p{t) = U{t)poW{t), we get a matrix representation £{t) = U*{t) (g) U{t) 
from the identity vec(Ai?C) = {C^ ^ A)vec{B) for any matrices A,B,C [52]. Here * denotes complex conjugation 
and ^ denotes transposition. The Schrodinger (von Neumann) equation describing evolution by a Hamiltonian H{t) 
translates to 



dt 



£{t) 



i(g> H{t) ^ H^it)(g)i £{t) 



':H{t)£{t), 



(Dl) 



with initial condition £{0) — I. In this appendix, we develop the generalized DHEOM that describs the reduced 
equations of motion for these propagators. 

Similar to the density matrix derivation, we begin by defining £{t,fl) = £ii{t)P{Q,t) where £(i is the conditioned 
propagator £nit) = £{t\fl{t)). This quantity satisfies a propagator version of the stochastic Liouville equation: 



d£{t,n) 
Ft 



£it,n) + rn£it,n), 



(D2) 



where we have split the Hamiltonian into a time- independent component and a Markov fiuctuating component: 
v. = Ho + ^{t)V . As with the density matrix, the averaged propagator is given by the marginal: 



£{t) 



dfi £{t,n). 



(D3) 



When Q{t) is a diffusive process with a generator possessing polynomial eigenfunctions, we expand P(ri, t) in generator 
eigenfunctions /„(fi) as 



£{t, n) = £n{t)P{n, t)^Yl ^n{t)U^). 



(D4) 



n=0 
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Then utilizing the three-term recurrence relations that the orthogonal polynomial eigenfunctions follow, we arrive at 
the following classical hierarchical equations of motion for the propagator: 





i 

n 




Bo 
Ao 


|4W = - 


i 

n 




Bn 
An 



V +Ao 



£o{t) 



i Ci 

hAi 



V 



£n{t)-i^V£nMt)- 



hA. 



n+l 



hAn- 



-V£n-lit), 



n> 0, 



(D5) 



where A„, An, Bn and C„ are eigenvalue and eigenfunction recurrence relation coefBcients defined identically as in the 
main text. The terminator equation for this DHEOM is explicitly: 



dt 



£N{t) 



Ho 



V + A. 



fiv(t) 



1 



C 



N+l 



fl^ Ajv+i^Ar+iA^r 



VV £N{t) 



1 



HA 



N-l 



V £N-i{t). 



(D6) 



